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Based on the nonrelativistic QCD factorization approach, 0{a a v 2 ) corrections to J/ip plus r) c 
production in e + e~ annihilation at yfs = 10.6 GeV is calculated in this work. The numerical results 
show that the correction at a s v 2 order is only about a few percent for the total theoretical result. 
It indicates that the perturbative expansions for the theoretical prediction become convergence and 
higher order correction will be smaller. The uncertainties from the long-distance matrix elements, 
renormalization scale and the measurement in experiment are also discussed. Our result is in 
agreement with previous result in ref [1]. 
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Study on heavy quarkonium decay and production is a very important and interesting issue to understand quantum 
chromodynamics (QCD), the fundamental theory of strong interactions. Many experimental and theoretical reserches 
have been performed since the discovery of the J/ip charmonium meson in 1974 followed by the T bottomonium 
meson in 1977, for reviews see Ref. [2]. In experimental side, it is easy to detect J/ip and T signal. In theoretical side, 
Qh> quarkonium bound states offer a solid ground to probe QCD, due to the high scale provided by the large mass of the 
heavy quarks, which make the QCD factorization possible in the related calculation. To explain the large discrepancy 
on the transverse momentum distribution of chromonium hadroproduction between the experimental measurement 
and theoretical prediction as well as to arrange the infrared divergence cancellation in p-wave quarkonium related 
calculation, the nonrelativistic QCD (NRQCD) factorization approach[3] has been introduced. It allows consistent 
theoretical prediction to be made and to be improved perturbatively in the QCD coupling constant a s and the 
| heavy-quark relative velocity v in heavy quarkonium rest frame. 
\Q . In last five years, most of the important theoretical studies on heavy quarkonium based on NRQCD are calculated 
f*"*. ' at next-to-leading order (NLO) of QCD and many of them are also calculated at next-to-leading order of v. Among 
them, the J/ip polarization puzzle at hadron colliders is still unclear after the important progresses at QCD NLO [4], 
It seems that the inclusive J/ip production at B-factories can be explained by just color singlet contribution at QCD 
t-H [ NLO [5, 6], but it causes the problem for the color-octet long distance matrix elements [7]. The theoretical calculation 
with NLO QCD and relativistic correction can cover the experimental measurements on exclusive double chromonium 
production at B-factories although the corrections are very large. For theoretical prediction based on perturbative 
expansion, the convergence of the expansion is a very important issue. Therefore, it is important to test the calculation 
' at higher order when the NLO correction is large. Usually, higher order calculation is much more complicate, so far 
, there are only a few simple processes whose 0(a s v 2 ) corrections are calculated [8-10]. 

For the exclusive double chromonium production at B-factories, its higher order calculation is studied, so we 
give a detailed review on it. The exclusive production cross section of double charmonium in e + e~ — > J/ipr] c at 
v/s = 10.6 GeV measured by Belle [11, 12] is a[J/ip + tj c ] x BV°[> 2] = (25.6 ± 2.8 ± 3.4) fb and by BABAR [13] 
is a[J/ip + rj c ] x B ,la [> 2] — (17.6 ± 2.8+3, i*) fb, where B Vc [> 2] denotes the branching fraction for the rj c decaying 
into at least two charged tracks. Meanwhile, the NRQCD LO theoretical predictions in the QCD coupling constant 
a s and the charm-quark relative velocity v, given by Braaten and Lee [14], Liu, He and Chao [15], and Hagiwara, 
Kou and Qiao [16] are about 2.3 ~ 5.5 fb, which is an order of magnitude smaller than the experimental results. 
Such a large discrepancy between experimental results and theoretical predictions brings a challenge to the current 
understanding of charmonium production based on NRQCD. Many studies have been performed in order to resolve 
the problem. From treatments beyond NRQCD, Ma and Si [17] treated the process by using light-cone method, a 
similar treatment was performed by Bondar and Chernyad [18] and Bodwin, Kang and Lee [19], possible contribution 
from intermediate meson rescatterings was considered by Zhang, Zhao, and Qiao [20], it was also studied in the 
Bethe-Salpeter formalism by Guo, Ke, Li, and Wu in Ref [21]. Based on NRQCD, Braaten and Lee [14] have shown 
that the relativistic corrections would increase the cross section by a factor of about 2, and the NLO QCD correction 
of the process has been studied by Zhang, Gao and Chao [22] and Gong and Wang [23], which can enhance the cross 
section with a K factor (the ratio of NLO to LO) of about 2, again the relativistic corrections have been studied by 
Bodwin, Kang, Kim, Lee and Yu [24] and by He, Fan and Chao [5], which is significant. More detailed treatment, 
such as including the resummation of a class of relativistic correction, has been taken into consideration by Bodwin 
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and Lee and Yu [25]. In another way, Bodwin, Lee and Braaten [26] showed that the cross section for the process 
e + e~ — > J/tp + J/tp may be larger than that for J /tp + r\ c by a factor of 1.8, in spite of a suppression factor a 2 /a 2 that 
is associated with the QED and QCD coupling constants. They suggested that a significant part of the discrepancy 
of J/tp + rj c production may be explained by this process. Hagiwara, Kou and Qiao [16] also calculated and discussed 
this process. In 2004, a new analysis of double charmonium production in e + e~ annihilation was performed by Belle 
[27] based on a 3 times larger data set and no evidence for the process e + e~ — >■ J/tp + J/tp was found. Both the 
NLO QCD corrections and relativistic corrections to e+e~ — > J/tp + r\ c give a large K factor of about 2. It is obvious 
that these two types of corrections to e+e~ — > J/tp + J/tp should be studied to explain the experimental results. 
In fact, they have been studied by Bodwin, Lee and Braaten for the dominant photon-fragmentation contribution 
diagrams [28]. The results show that the cross section is decreased by K factor of 0.39 and 0.78 for the NLO QCD and 
relativistic corrections respectively. A more reliable estimate, 1.69 ± 0.35 fb, was given by Bodwin, Lee, Braaten and 
Yu in ref. [29]. And light-cone method is used in ref. [30] by V.V. Braguta. Gong and Wang performed a complete 
NLO QCD calculation on e + e~ — > J/tp + J/tp [31] and the results show that the cross section would be much smaller 
than the rough estimate in Ref. [28]. Therefore it is easy to understand why there was no evidence for the process 
e + e~ — > J/tp + J/tp at B-factories. 

It is easy to sec that both the QCD correction (a s ) and relativistic correction (v 2 ) are very large for e+e~ — > J/tpr] c 
at B-factory energy, and with these corrections the experimental measurement can be explained. Therefore it is natural 
to ask the question, how is the situation for the higher order corrections beyond a s and v 2 correction ? a 2 correction 
is very difficult to do, but recent progress make it available to do a s v 2 correction already. It is very interesting to see 
that the a s v 2 correction, given in a recent work [1], is a small contribution. It convinces us in some sense (with a 2 
correction absent) that the double expansions in NRQCD converges quite well on this problem. Since the calculation 
is quite complicate and plays an important role to convince us the convergence on the theoretical predication which 
can explain the experimental data, in this paper we performed an independent calculation on it by using the our 
package Feynman Diagram Calculation (FDC) [32] with the built-in method to calculate relativistic correction. The 
remainder of this paper is organized as follows. Base on the NRQCD frame, we briefly introduce theoretical formulism 
for the calculation of heavy quarkonium production and give the corresponding results in perturbative NRQCD in 
Sec. II. The details in perturbative QCD are summarized in Sec. III. We give the numerical results of a s v 2 corrections 
and some discussion in Sec. IV. Finally, in Sec. V, we present a brief summary. 

II. NRQCD FACTORIZATION FORMULA UP TO v 2 ORDER 

According to NRQCD effective theory, the production of the charmonium are factorized into two parts, the short- 
distance part and the long-distance part. The long-distance parts are related to the four fcrmion operators, character- 
ized by the velocity v of the charm quark in the meson rest frame. The long-distance matrix elements can be estimated 
by lattice calculations or phenomenological models, or determined by fitting experimental data. The production cross 
section up to v 2 order is expressed as 

<j(e + e~ -> J/tp + r] c ) = (c 00 + c w {v 2 ) J/il + c i(w 2 ) r ,J(Oi) r , c (d) J/v , (1) 
with the long-distance matrix elements being defined by using related operators as 

tf)jH> = T/nT , <Oi>j/* = <0|xW(at a^VVxlO), (2) 

for J/tp and 

{y2) ^ = Sin\ ■ = (0|xV(4c^c)^x|0), (3) 

m Z{ U l)r,c 

<Pi>,e = (0|i[x t V'«a T ,JV t (-^) 2 X + X t (-^) 2 ^««,JV' t x]|0), 

for ?7 C where m c is the charm quark mass. It is the basic point that the NRQCD factorization for hadron related 
process will also hold when the hadron state are replaced by QQ states with exactly the same quantum numbers 
as the corresponding hadron state. In this way, the short-distance coefficients c o,c i and c 10 can be obtained in 
perturbative calculation through the matching condition, and they are calculated up to QCD next-to-leading (NLO) 
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order. In order to obtain the short-distant coefficients, the matrix elements of the operators for quantum states need 
to be calculated perturbatively, and there are 

(0 1 ) 1Sa = 2N c (2E qi ) 2 , (O^s, = W c (2E q2 ) 2 (4) 



where there are N c = 3 for S£/(3) group and E q = yj m 2 + q 2 . From the NRQCD effective Lagrangian, we could 
easily get the Feynman rules. Therefore we have calculated order a s v 2 corrections to the leading order (Oi}2s+i Ss 
in perturbative NRQCD with the dimensional regularization and defined the renormalization constants Z^ s of the 
operator by using the MS scheme [3, 33]. 

r 7 ls 4:a s C F [i 2 1 q 2 

SZ = -__(_)(_+ In 47t-7e)-2 (5) 
(Oi>5 + i S . = [1 + ^^(4) e (- f + ^tt - 7e)4(Oi) 2s+1 s s . (6) 

07T r^A c 



(V 1 )2 S+1Ss =q 2 (0 1 ) 2 , +15s . (7) 
At last we could easily give the perturbative NRQCD results. 

.(e+e- - QQeSl) + QO^o 1 ))^^ = {coo + ^[c 10 + )^ + h*r - 7E )e° ] + ^ 

[ CQ1 + ^£(^ r( I +ln4 7 r- 7E )4]}192(iV c ^ 1 ^ 2 ) 2 (8) 
3tt p A e 



III. DETAILS OF PERTURBATIVE QCD CALCULATION 



For a Q(p)Q(p) quantum state, we denote P as the total momentum and q as the relative momentum between Q 
and Q pair. Therefore, there are 

P=\p + q, P=\p-q- (9) 
p 2 =p 2 =m 2 Q , P 2 =AE 2 , E q = ^m 2 Q + q 2 

where tuq is the mass of the heavy quark Q, and the Q and Q are on their mass shells. 

To do the perturbative calculation in related process for the quantum states, we could obtain the projectors for each 
quantum states. The spin-singlet and spin-triplet components of each QQ state can be projected out by making use of 
the spin projectors. After multiplying corresponding Clebsch-Gordan coefficients to the spin component of the outer 
product of the spinors for each QQ pair, one can find that LTi and II3 are the spin-singlet and spin-triplet projectors 
of the QQ production, respectively. The spin projectors that are valid to all orders in the relative momentum can be 
found in Rcfs[34]. 



Ill 

n 3 



1 



AV2E(E + mq) 
1 

4V2E(E + m Q ) 



(p-m Q ) l5 (P+2E)(p' + m Q ), 
(P-m Q )f(\)(r+2E)(p' + m Q ), 



(10a) 
(10b) 



where II 1 and II3 are projectors for spin and spin 1 s-wave quantum states respectively, and e*(A) is the polarization 
vector of the spin-triplet state. 

For process e+ (pi ) e~ (p 2 ) -> Q ( f - qi ) Q ( f + qi ) ( 3 Sj ) + Q ( f - - q 2 ) Q ( ^ - q 2 ) ( 1 So ) > thc production matrix element 
is expressed as 



M (e+e- -> QQ( 3 Sl) + QQ( 1 S 1 )) = e„(S z )A»( qi ,q 2 ) 



e^(S z ) A" 



+ q 2 iap d 3 A^ 
91=0,92=0 6 dq^dqf 



+ Si I<*0 

9l =0,92=0 6 



d 3 A M 



dq^dq^ 



91=0,92=0 



+ 0(qlq 4 2 ) 



(11) 
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where we have used the following relation 

As for the expansion of q, we should consider the effect that the external momentum and polarization vector may be 
the implicit function of q. From the momentum conservation and on-shell conditions, p\ = 4£? 2 i ,p| = 4£? 2 2 , we could 
find that p 3 , p^ are implicit functions of qi, q2 respectively. However, it is obvious that the short-distance coefficients, 
to be obtained in the perturbative calculation, are functions of the independent variables which are the invariant mass 
s of the e+ and er system and cos 9. 9 is the angle between J/ip and the electron. Where s and cos9 are independent 
of the relative momentum q. 

Since the final results are Lorentz invariance and irrelevant to the reference frame, we choose to do the calculation 
in the center-of-mass of this system where p\ +P2 = P3 +Pa = (y/s, 0, 0, 0) is the explicit expression of the momentum 
conservation. Therefore the following results are obtained: 



dp 3 „ dp 4 dp 3 dp 4 dp 3 

X~2"--P3 = 2, T - s -p 4 = 0, + =o, _. p4 = 0. (13 

dq( dq( dq( dqf dq( 

We choose two vectors n = (0,ri) and r 2 = (0, ¥2) with rt and r| being unit vectors, while rt, r| and p% are 

d£3 
dq? 

dp 3 
dq\ 

dp 3 n dp 3 

dqf dqf 

together with previous conditions in Eq.(13), we can easily obtain the solution 



perpendicular to each other, i.e n • r 2 = 0,p 3 ■ n = 0,p 3 • r2 = 0. Then vector — ^ can be expressed as linear 

if 

combination of four independent vectors as ^| = aip 3 + a 2 p4 + a 3 ri + a^r 2 - From the following conditions 

2-ri=0, ^T2=0 (14) 



dp 3 -2p| 2p 3 • p 4 



+ 71 — TT2 ( 15 ) 



dq? (P3 • P4 ) 2 - p 2 pl (P3 • P4 ) 2 - P 2 p1 

^3 Q„ ^\ -^T-^fl-i Vicili pHir \ ■fViavti euro vnlof inn ^ ^. 

1? 



For the e*(A), the polarization four-vector of the \QQ( 3 Si)} with helicity A, there are the relation — — — = since 

dq x 



9 is independent of the relative momentum q. It is easy to obtain 

^qr" P3 -~dqf e (0) ' ^r' e(0) - ' ^qr' e(1) -°' -^T" 6 ^- - (16) 

Therefore, we obtain the relation between the polarization four-vector and q as 

dp3 * m 

de* (A) dq? ■ 6 lAJP3 _ -2p 3 ■ p iPA ■ e* (X)p 3 



dq? p\ {(p 3 ■ Pa) 2 - pIpDpI ' 

The treatment about q 2 is similar to these. 

Considering the two body phase space, we need to expand it. 

dT= I <l c,»<)- 2 



(17) 



-> A 1 /2( s ,4E? l ,4£7?,) 222 ~t 

where \p \ = -£= — , X(x, y,z) — x + y + z — 2(xy + yz + xz). Only \p | need to be expanded since 

2yfS 

cos 9 and s are independent of q\ and g 2 • Then there is 



where = — — -. We could square the amplitude, integrate over the phase space, and expand in powers 

2^s 

of q in order to obtain the desired perturbative result 
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a(e+e- -+ QQ( 3 Sl) + QQfSfo) 



pcrtQCD 



M 



(18) 



Most of the steps in this section are realized in a small program in FDC package, and the final Fortran source for 
numerical calculation are prepared by using FDC package together with the small program for q 2 expansion. 

Since there is no 0(a s ) real process in NLO, we only need to calculate virtual corrections. Dimensional regularization 
has been adopted for isolating the ultraviolet (UV) and infrared(IR) singularities. UV divergences are cancelled upon 
the renormalization of the QCD gauge coupling constant, the charm quark mass and field, and the gluon field. A 
similar renormalization scheme is chosen as in ref. [35] except that both light quarks and charm quark are included 
in the quark loop to obtain the renormalization constants. The renormalization constants of the charm quark mass 
Z m and field Z 2 , and the gluon field Z 3 are defined in the on-mass-shell(OS) scheme while that of the QCD gauge 
coupling Z g is defined in the modified-minimal-subtraction(MS) scheme: 



sz: 



OS 



SZ os 



2 — 



6Z°* 



5Z™ S = - 



-3C F 



Us 
An 



1 , Anfj 2 
7E + In =- 

euv mi 



1 



A 
3 



•0(e) 



2 Attn 2 

3 7B + 31n— tL+4 + 0(e) 

em mi 



a s 
An 



($-2CU)(- -) -tT F [ — - 7 £ + ln^ 

\£uv e IR J 3 \e uv mj. 



2 4tt 



1 



(_L__L)_W_L 

j E + 1h(4tt) + 0(e) 



+ 0(e) 



(19) 



where 7# is Euler's constant, /3o = TfCU — ^Tprif is the one-loop coefficient of the QCD beta function and rif is 
the number of active quark flavors. There are three massless light quarks u,d,s, and one heavy quark c, so nf=A. 



In SU(3) C , color factors are given by T F = \,C F = |,C A = 3. And /3' = /3 + (4/3)T F = (11/3)CU - (4/3)7>n z/ 



2,<^F — 3, — ^liiu 

nif = nf — 1 = 3 is the number of light quarks flavors. Actually in the NLO total amplitude level, the terms 

OS 



where 

proportion to SZ3 



cancel each other, thus the result is independent of renormalization scheme of the gluon field. 



IV. RESULTS 



The final results are obtained by using the matching method with the UV and IR divergences being cancelled. 

°" = &LO + °NLO(u s ) + &NLO(v 2 ) + &NLO(av 2 ) (20) 

&LO,VNLO(a s ),< J NLO(v 2 ),o'NLO(av 2 ) arc the contributions from the leading order, the next leading order in a s , the 
next leading in v 2 and the next leading in av 2 . Then the production rate up to 0(a s v 2 ) order is expressed as 



a = 



_ 81927r J Cg.e>;( t t r )q^(l-4r) 



S/2 



9N 2 s 



■^^[/^oln^/iW + f ln£ + / 4 (r)] + ^< [/3 In £-/ a (r) + f In £ 



+■ Mr)] 
Mr)]} 



(21) 



where there are e c 



and u r is the renormalization scale. 



4 r — 4m c f _ 9-74r+80r 2 f ( \ _ ll-82r+80r 2 

3, r - — , fi{r) - 6 (i-4r) ' - — 6 (i-4r) — 

Therefore, the obtained analytic expressions of the v 2 correction is in agreement with that in the paper [1] . At the 
same time, the analytic expression of / 3 (r) in the results of the a s correction is also in agreement with that in the 
paper [23]. Since the analytic expressions of fi{r) and /s(r) in the 0(a s v 2 ) correction are so lengthy that we just 
give the numerical results for them. In the numerical calculation, there are 



4x1. 4 2 



12.358, U = 3-8382, f 5 = 3.2537, for r = TW , 
11.806, ,U = 2.0543, f 5 = 2.6668, for r 



4x1.5-' 
10. 58 2 



h = 0.97466, f 2 = 1.3080, f 3 
h = 0.87465, h = 1-2098, / 3 

And we take y/s = 10.58 GeV and fj,\ = m c . The running strong coupling constant is evaluated by using the 
two- loop formula with A^ = 0.338 GeV as used in Ref [23]. Our results are presented in the table I with parameters 
given in table caption. The results are in agreement with that in ref [1]. The contribution from the 0(a s v 2 ) order 
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a s (fi r ) 




°"jVLO(a s ) 


a NLO(v 2 ) 


a NLO(av 2 ) 


cr 


a,(a£) = 0.211 


4.381 


5.196 


1.714 


0.731 


12.273 


a s (2m c ) = 0.267 


7.0145 


7.367 


2.745 


0.245 


17.372 



TABLE I: With the follow parameters: a( v / s) = 1/130.9, (Oi)j/^ = 1.161GeV 3 , (d)^ = 0.387GeV 3 , m c = 1.4 GeV, 
{ v2 )j/i> = 0.223 , {v 2 ) Vc = 0.133, We give the cross sections with different renormalization scale fx. Their units are fb. 




K (GeV) 



FIG. 1: The cross section as a function of the fj, r at y/s = 10.58 GeV. The black and blue solid curves are the cross sections 
in the m c = 1.4 and m c = 1.5 respectively. The red and green bands represent the measured cross sections by the Belle and 
BaBar experiments, with respective systematic and statistical errors. 



is small. There are some differences for the results in the table II if the long-distance matrices and QED coupling 
constant are chosen as in Ref [5], the correction at 0(a s v 2 ) order is also small, we also give the relation of the cross 
sections and the fx r in the FIG.l. There are about 10 percents differences in the total cross sections between m c = 1.5 
and m c = 1.4. We find the uncertainty of the total cross sections from m c is not small. If we choose [i r = 2m c , we 
could give the relations of the ratios of different parts to total cross sections with the yfs in the FIG. 2. We will find 
that the contributions from 0(a s ) and 0{a s v 2 ) become important and the one from LO becomes small when y/s is 
large, but there are just about 10 percents contribution in 0(a s v 2 ). The contribution from the 0(a s v 2 ) order is small 
once again. 

V. SUMMARY 

In this work we have calculated the 0(a s v 2 ) correction in detail for the processes e + e~ — > J/ip + r\ c within the 
frame of NRQCD. The result at 0(a s v 2 ) order give about 6 percent contribution to the total theoretical prediction 
while the 0(a s ) correction and 0(v 2 ) are about 40 percent and 14 percent contribution respectively. It indicates 
that the convergence in the double perturbative expansions in QCD a s and relativistic v 2 are very well for the 
theoretical calculation on the production rate of the process e + e~ — > J/ip + rj c . Up to 0(a s v 2 ) order, the theoretical 



m 


a s (fj, r ) 




® N LO{a s ) 


a NLO(v 2 ) 


a NLO(av 2 ) 


a 


1.5 


= 0.211 


5.973 


6.645 


1.335 


0.416 


14.369 


1.5 


a s (2m) = 0.259 


9.000 


8.771 


2.011 


-0.017 


19.726 


1.4 


«.(f) = 0.211 


6.526 


7.754 


1.591 


0.667 


16.538 


1.4 


a,(2m c ) = 0.267 


10.450 


10.989 


2.548 


0.1989 


24.185 



TABLE II: In the follow parameters: a(y/s) = 1/137, {Oi)j/^ = 1.719GeV 3 , {Oi) Vc = 0.432, GeV 3 , {v 2 }.;/^ = 0.090 , {v 2 ) Vc = 
0.119,, We give the cross sections with different m and renormalization scale /i. 
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prediction with quite large uncertainty from charm quark mass and renormalization scale can describe the experimental 
measurement. 

This work is supported by the National Natural Science Foundation of China (Nos. 10979056 and 10935012), in 
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